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This  work  shows  that  optimal  compact  storage  of 
coefficient  matrices  affords  a  significant  reduction 
in  core  storage  requirements  over  banded  storage 
schemes.  The  resulting  savings  enables  in  core  fi- 
nite element  solutions  of  large  systems  not  otherwise 
possible.  It  is  shown  that  Gears  method  for  the  stiff 
system  of  a  nonlinear  reactor  dynamics  problem  is  not 
as  efficient  as  Crank-Nicolson  integration  because  of 
substantially  greater  core  requirement,  despite  its 
superior  tracking  ability.  A  remedy  in  the  form  of  a 
modified  implicit  version  of  Gears  method  with  a  sig- 
nificant reduction  in  core  requirement  is  shown  to  pro- 
vide the  same  excellent  accuracy  as  Gears  method.  Com- 
parisons between  the  modified  Gear  method  and  the 
Crank-Nicolson  method  show  the  relative  advantages 
and  disadvantages  of  each.  Finally,  it  is  shown 
that  although  the  nonlinearity  encountered  in  this 
problem  can  be  treated  directly,  a  linear  approxima- 
tion of  the  nonlinear  term  affords  a  substantial 
reduction  in  core  requirement  with  a  relatively  small 
cost  in  accuracy. 
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AN  OPTIMAL  COMPACT  STORAGE  SCHEME  FOR 
NONLINEAR  REACTOR  PROBLEMS  BY  FEM 

by 

D.  Salinas,  D.  Nguyen  and  R.  Franke 

Abstract  -  This  work  shows  that  optimal  compact  storage  of  coefficient 
matrices  affords  a  significant  reduction  in  core  storage  requirements  over 
banded  storage  schemes.  The  resulting  savings  enables  in-core  finite 
element  solutions  of  large  systems  not  otherwise  possible.  It  is  shown 
that  Gears  method  for  the  stiff  system  of  a  nonlinear  reactor  dynamics 
problem  is  not  as  efficient  as  Crank-Nicolson  integration  because  of 
substantially  greater  core  requirement,  despite  its  superior  tracking 
ability.  A  remedy  in  the  form  of  a  modified  implicit  version  of  Gears 
method  with  a  significant  reduction  in  core  requirements  is  shown  to  pro- 
vide the  same  excellent  accuracy  as  Gears  method.  Comparisons  between  the 
modified  Gear  method  and  the  Crank-Nicolson  method  show  the  relative 
advantages  and  disadvantages  of  each.  Finally,  it  is  shown  that  although 
the  nonlinearity  encountered  in  this  problem  can  be  treated  directly,  a 
linear  approximation  of  the  nonlinear  term  affords  a  substantial  reduc- 
tion in  core  requirement  with  a  relatively  small  cost  in  accuracy. 

Statement  of  the  Problem 

The  finite  element  method  (FEM)  has  been  successfully  applied  to  the 
problem  of  a  space- time  nonlinear  reactor  with  temperature  dependent  feed- 
back (1_,2J.  However,  the  effective  numerical  treatment  of  the  resulting 
system  of  stiff  O.D.E.  called  for  further  investigation.  This  work  com- 
pares alternative  computational  schemes  in  a  FEM  solution  of  a  nonlinear 


nuclear  reactor  problem. 

A  test  case  is  chosen  which  consists  of  a  fully-reflected  fast 
reactor,  subject  to  super  prompt  critical  conditions.  By  neglecting  the 
delayed  neutrons  and  using  the  prompt  feedback  model,  the  following  non- 
linear reactor  dynamics  equation  is  obtained  in  (r,z)  space  in  the  core, 

1  ||  (r,z,t)  _  Da2^  +  xz^   .  w^2  =  0  (1) 

A  similar  equation  without  the  nonlinear  feedback  term  results  in  the 

reflector.     A  detailed  development  of  these  equations  is  given  in  ( 1_) . 

In  accordance  with  Galerkin's  method,  Eq.   (1)   is  transformed  into 
the  implicit  system  of  N  ordinary  differential  equations, 

e    (A^ijttJ  -  B^j)  +    z        i      CijkVk  =  0  (2) 

j=l  j=l     k=l 

i,  =  1,  ...,  N 
together  with  appropriate  initial  conditions,  ^.(0),  i  =  1,  ...,  N.  The 
computational  efficiency  of  any  finite  element  program  depends  on  the  al- 
ternative schemes  available  for  selection  by  the  investigator.  The  pre- 
sent work  considers  the  relative  merits  of  various  integration  methods, 
algebraic  equation  solvers,  array  storage  methods,  as  well  as  the  effects 
of  linearization. 

Linearization 

A  well  known  procedure  for  treating  nonlinear  terms  in  FEM  is  by  a 
"linearization"  of  the  nonlinear  terms.  Linearization  here  refers  to  the 
process  whereby  a  nonlinear  term  F  [if>-(t)]  is  approximated  by  F  [^(t-At)], 
where  ^(t-AtJis  the  known  value  of  ^  at  the  previous  time.  For  the  par- 


N  N 

icular  quadratic  nonlinear  term  z  Z       C  ..  ^.(t)^.  (t)  of  Eq.  (2),  a 

j=l  k=l   1JK  J    K 
better  linear  approximation  is 

N 

z  C*.  \l>At),     where 
j=l  1J  J 


Cij  =  &     C1jk*k(t-"At)  (3) 


It  should  be  noted  that  although  the  nonlinear  quadratic  term  eeC...^ -ij;. 

I J  K  J  K 

can  be  treated  directly,  and  indeed  it  was  in  reference  (1_) ,  the  linear 
approximation,  Eq.  (3),  results  in  a  substantial  reduction  in  both  array 
storage  and  CPU  time  for  a  relatively  small  expense  in  accuracy,  less  than 
5  percent  for  the  cases  run.  The  savings  in  storage  and  CPU  time  depends 
on  the  integration  method  and  the  algebraic  equation  solver  employed,  as 
well  as  the  method  of  array  storage. 

By  far  the  most  common  (but  not  optimal)  method  of  array  storage  is 
the  band  storage  scheme,  whereby  only  the  bandwidth  (or  half  bandwidth  in 
the  case  of  symmetric  matrices)  terms  in  a  matrix  are  stored.  Since  ma- 
trices obtained  in  FEM  analyses  are  always  banded,  banded  storage  is  always 
advantageous.  Table  1  lists  the  storage  requirements  and  the  storage 
savings  accomplished  by  linearization  of  symmetric  and  nonsymmetric  C-.  . 
Coefficient  a  is  the  number  of  bytes  per  word  for  the  particular  computer 
being  utilized.  For  an  IBM360  computer,  a=4  for  single  precision,  and  8 
for  double  precision.  N  is  the  number  of  unknowns  and  n  is  the  total 
bandwidth. 


TABLE  1 
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'ijk 
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Reduction 


non  sym.,non  banded    aNN 


aN 


aN2(N-l) 


non  sym. , banded 


aNn' 


aNn 


aNn(n-l ) 
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sym., non  banded    — ^<N+1 ) (-^  +  — g — )   — 2 — 


g(N+l)N  /2N+1   1* 
9     \    c    ~  T) 


sym. , banded 


;aj(W)^+-a^j   sN(n+1) 


_{n+l)(__.  _) 


It  should  be  pointed  out  that  some  additional  reduction  can  be  obtained 
for  symmetric  matrices  by  discarding  the  n  g  '     zeroes  in  the  lower  right 
triangle  of  the  banded  N  x  n  Ct.  matrix.  In  most  cases  n  is  small  relative 
to  N  and  the   g  '  additional  reduction  is  not  worth  the  effort. 

Table  2  gives  an  idea  of  the  savings  achieved  using  banded  storage 
for  some  typical  regular  rectangular  grids  as  shown  in  Figure  1.  The 
results  are  associated  with  the  symmetric  C  ..  and  £f.  matrices  of  the  non- 
linear  reactor  dynamics  problem  considered  in  this  work,  in  which  linear 
triangular  elements  were  utilized  with  single  precision  on  an  IBM360  com- 
puter. (R-l)  is  the  number  of  rows  and  (S-l)  the  number  of  columns  in  the 
rectangular  grid,  with  S  >  R. 


TABLE  2. 


j  Grid" 

Cijk 

c*. 

Reduction 

-  R=5  S=8 

;  N=40  n=3 

i 

5,600 

1,120 

4,480 

!  R=10  S=20 
j  N=200  n=23 

1 

80,000 

9,600 

70,400 

'  R=20  S=40 
N=800  n=43 

1 

,056,000 

70,400 

985,600 

R=40   S=60 
;N=2400  n=83 

11 

,424,000 

403,200 

11 

,020,800 

Optimum  Compact  Storage  of  Coefficient  Matrices 

Since  the  interpolation  function  for  the  i   unknown  is  non-zero  only 
over  those  elements  which  contain  i  (i.e.,  the  basic  functions  in  FEM 
have  local  support  only),  the  matrices  resulting  from  FEM  are  both  banded 
and  sparse.  The  sparce  nature  arises  from  the  fact  that  the  unknowns 
surrounding  the  i   unknown  are  not  consecutively  numbered.  Common  prac- 
tice takes  advantage  of  the  banded  nature  of  these  matrices  via  a  banded 
storage  scheme,  which  stores  an  NxN  matrix  as  an  Nxn  matrix  where  N  is 
the  number  of  unknowns  and  n  is  the  total  bandwidth  of  the  matrix. 
Using  banded  storage,  special  care  in  the  numbering  of  unknowns  can  affect 
a  substantial  decrease  in  n.  As  a  simple  example  of  this,  consider  the 
regular  rectangular  grid  of  Figure  1,  having  (R-l)  rows  of  elements,  and 
(S-l)  columns  of  elements.  For  linear  triangular  elements,  sequential 
vertical  numbering  results  in  n  =  2R  +  3,  whereas  sequential  horizontal 


numbering  gives  n  =  2S  +  3.  The  saving  achieved  by  vertical  numbering, 
assuming  R  <  S,  is  2N(S-R)  words,  or  2aN(S-R)  bytes,  where  a  is  the  num- 
ber of  bytes  per  word.  For  symmetric  matrices  the  previous  numbers  are 
halved. 

These  results  may  be  generalized  for  triangular  elements  of  order  t, 
For  the  rectangular  grid  of  Figure  1,  the  bandwidth  for  triangular  ele- 
ments  with  t   order  polynomial  interpolation  is 


n  =  6Rt  -  4R  -  2t  +  5 


for  non  symmetric  matrices,  and  the  half  bandwidth  is 


ns  =  3Rt  -  2R  -  t  +  3 


(4) 


(5) 


for  symmetric  matrices.  Equations  (4)  and  (5)  are  derived  on  the  basis  of 
condensing  out  the  interior  element  nodes.  Table  3  presents  the  evalua- 
tions of  these  formulas  for  linear,  quadratic  and  cubic  elements  for  the 
rectangular  grid  of  Figure  1. 

TABLE  3. 


_____ 

t 

n 

ns 

1 

2R  +  3 

R  +  2 

2 

8R  +  1                     4R  +  1 

3 

14R  -   1                      7R 

For  banded  systems  of  linear  algebraic  equations,  direct  methods  of 
the  Gaussian  elimination  type  are  most  often  utilized,  being  preferred  over 
iterative  methods  such  as  successive  over-relaxation  (SOR).  Depending  on 
the  particular  application,  there  are  points  in  favor  of  each  method.  In 
the  case  of  linear  equilibrium  (elliptic)  boundary  value  problems  Gaussian 
elimination  has  the  following  advantages: 

1.  Direct  methods  allow  the  use  of  iterative  refinement. 

2.  If  more  than  one  right  hand  side  must  be  processed,  the  initial 
cost  of  decomposition  must  be  expended  only  for  the  first  solu- 
tion. 

For  transient  (parabolic  or  hyperbolic)  initial-boundary  value  problems, 
or  multi-step  nonlinear  equilibrium  problems,  SOR  has  the  following 
advantages: 

1.  The  large  number  of  iterations  for  a  solution  at  the  first  step 
is  expended  only  once,  succeeding  solutions  require  far  fewer 
iterations  since  a  good  estimate  of  the  solution  is  then  avail- 
able. 

2.  No  working  space  is  required  by  SOR  and  therefore  matrices  may 
be  compacted  to  maximum  density. 

It  is  this  last  point  which  suggests  the  replacement  of  the  banded  storage 
schemes  by  optimum  compact  storage  (OCS)  which  will  be  described  shortly. 
The  argument  in  defense  of  direct  methods  that  storage  is  becoming  in- 
creasingly abundant  is  a  spurious  one  which  impedes  investigation  of 
larger  systems,  especially  for  three  dimensional  problems.  It  will  be 
shown  here  that  optimum  compact  storage  of  coefficient  matrices  leads  to 
a  substantial  saving  of  core  for  large  systems.  A  description  of  OCS 


follows. 

Accounts  of  OCS  are  given  by  Tewarson  {3)   and  Gustavson  (4J .  The 
underlying  idea  behind  OCS  is  simply  to  store  only  the  non  zero  terms  of 
a  coefficient  matrix.  The  reduction  in  storage  requirement  achieved 
through  OCS  as  compared  to  banded  storage  is  illustrated  by  the  regular 
rect  angular  grid  shown  in  Figure  1.  It  is  striking  that  regardless  of 
R  and  S,  the  band  of  n  terms  surrounding  any  nodal  point  contains,  at  most, 
7  non  zero  coefficients  in  the  case  of  linear  triangular  elements.  For 
the  cases  of  quadratic  and  cubic  elements  the  maximum  number  of  non  zero 
coefficients  for  any  nodal  point  is  19  and  30  respectively.  Hence  the 
maximum  storage  requirements  for  the  regular  grid  of  Figure  i  for  linear, 
quadratic  and  cubic  triangular  elements  are  a(7N),  a(19N)  and  a(30N), 
respectively.  The  actual  reduction  is  not  simply  the  difference  between 
banded  core  requirement  and  OCS  core  requirement  for  two  reasons. 

1.  OCS  storage  can  be  achieved  with  less  than  a(yN)  where  y   is 
the  maximum  number  of 'non  zero  coefficients  and  a  is  the  number 
of  bytes  per  word. 

2.  Additional  vector  arrays  are  required  for  the  implementation  of 
OCS. 

Although  the  following  discussion  is  given  in  terms  of  the  linear 
NxN  C*.  •  array,  the  nonlinear  NxNxN  C  ..  array  is  treated  in  the  same  way. 

1  J  1  j  K 

Indeed  the  direct  treatment  of  the  nonlinearity  by  OCS  provided  the  most 
impressive  reduction  in  core  imaginable.  An  additional  word  on  this  will 
be  given  at  the  end  of  this  section.  It  is  well  to  keep  in  mind,  however, 
that  although  direct  treatment  of  the  nonlinearity  is  esthetically  pleas- 
ing, it  is  a  questionable  luxury  in  view  of  the  excellent  results  obtained 


by  linearization. 

The  implementation  of  OCS  requires  two  integer  array  vectors,  say 
ISTART  and  NAME,  and  a  vector  of  the  non  zero  coefficients  of  C,  which 
we  shall  call  CC.  The  i   integer  entry  in  the  (N+l)xl  ISTART  vector 

is  the  number  q.,  where  q.  =    Z     p.+l ,  and  p.  is  the  number  of  terms 

i=l 
th      .  J 

in  the  l   equation,  i.e.,  p.  is  the  number  of  nodes  surrounding  the 

+"  h  f*  h 

i   node.  ISTART,  then,  is  a  pointer  vector  whose  i   term  locates  the 

initial  position  in  the  CC  vector  of  the  contributing  coefficients  to 

the  i   equation.  The  Mxl  NAME  vector,  where  M  =  E  p.,  is  composed  of 

i=l 
N  successive  vector  blocks  of  variable  length  p.,  i=l,...,N.  The  p. 

integer  numbers  in  the  i   block  of  NAME  identify  the  p.  contributors  to 
the  i   equation.  The  Mxl  CC  vector  contains  the  real  non  zero  coeffi- 
cients of  the  NxN  matrix  C,  arranged  in  the  same  contiguous  block  arrange- 
ment as  the  NAME  vector.  The  jth  term  in  the  ith  block,  CC  (ISTART 
(I)  +  J-l)  is  coefficient  Cj.  where  K  =  NAME  (ISTART  (I)  +  J-l).  To  fix 
ideas  the  vectors  associated  with  the  simple  grid  shown  in  Figure  1  are 

ISTART  =  <  1   5   10   13   18   25   30   33   38   41  > 
NAME   =<  1   2  4  5 1  2  3  1  6  5 1 |  9  8  7  > 

CC    =  <  cn  ci2  C14  Cl 5 1  C22  C23  C21  C26  C25''  *  " '  C99  C98  C97> 
In  this  case  N=9  and  M=40. 

With  the  interior  nodes  of  an  element  condensed  out,  the  length 
(number  of  entries)  of  the  NAME  and  CC  vectors  for  the  regular  rectangular 

J.L. 

grid  of  Figure  1  with  triangular  elements  of  t   order  polynomial  interpo- 
lation, is: 


N 
Z 

i=l 


M  =  Z   p.  =  RS  (15t2  -  6t  -  2)  +  (R  +  S)  (-14t2  +  8t  +  2) 


+  (13t2  -  lOt  -  1) 


The  number  of  unknowns  is 


N  =  RS  (3t  -  2)  +  (R  +  S)  (2  -  2t)  +  (t  -  1) 


(5) 


(6) 


Table  4  gives  the  evaluation  of  Equations  (5)  and  (6)  for  linear,  quadratic 
and  cubic  triangular  elements  for  the  grid  of  Figure  1. 

TABLE  4. 


RS 


7RS  -  4  (R  +  S)  +  2 


4RS  -  2  (R  +  S)  +  1 


46RS  -  38  (R  +  S)  +  31 


7RS  -  4  (R  +  S)  +  2 


115RS  -  120  (R  +  S)  +  86 


Hence  the  core  requirements  for  banded  and  compact  storage  are: 


Np  =  aNn* 


n*  =  n  for  non  sym.  matrices 
n*  =  n/2  for  sym.  matrices 


N  =  aM  +  6  (M  +  N  +  1) 


(7) 


respectively,  where  a  and  3  are  the  number  of  bytes  per  word  for  real  and 

integer  numbers.  For  an  IBM360,  a  is  4  for  single  precision,  and  8  for 

15 
double  precision;  6  is  2  for  integers  less  than  2   and  4  for  integers 

15 
greater  than  2  .  Comparative  results  between  banded  and  OCS  storage, 


10 


using  the  formulas  from  Table  2  and  3,  for  single  precision  with  an  IBM360 
computer  for  symmetric  operators  on  some  typical  RxS  rectangular  grids  are 
given  in  Table  5. 

TABLE  5. 
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2,400 

9,401 

16,402 

i 
M    ! 

1,282 

8091 

19,486;    5,362 

I 

i 

34,551 

84,886 

16,402 

106,686 

264,086 

i 

ns    | 

i               1 

12 

41 

70j          22 

81 

140 

42 

161 

280 

> 

;nb  l 

9,600  121,524 

358,960:70,400 

998,244' 

3,002,720 

403,200 

6,054,244 

18.4x10° 

i 

Nc 

8,094:50,030 

119,482 

33,774 

j 
213,470 

520,042 

103,214 

658,920 

1,617,322 

NB"Nc 

1,506 

; 71 ,492 

239,478 

36,626 

784,774 

2,482,678 

299,986 

5,395,324 

16.8xl06 

X  Dif 

i 
1 

-15.7 

-58.5 

-66.7 

-52.0 

-78.6 

1 

-82.7 

-74.4 

-89.1 

-91.2 

The  notation  in  Table  5  is  as  follows: 

t  -  order  of  polynomial  interpolation 

N  -  number  of  unknowns 

M  -  number  of  terms  in  the  NAME  and  CC  vectors 

n  -  bandwidth  for  the  symmetric  matrix 


L  -  core  required  for  banded  storage,  single  precision  IBM360 
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N  -  core  required  for  OCS,  single  precision  IBM360 

Ng-N  -  reduction  in  core  storage  by  OCS,  single  precision  IBM360 

%   Diff  -  percent  difference  between  banded  storage  and  OCS 

In  Table  5,  there  is  no  intention  to  suggest  comparisons  for  a  par- 
ticular grid  with  respect  to  the  order  of  polynomial  interpolation.  Com- 
parisons with  respect  to  N,  the  number  of  unknowns,  are  more  meaningful. 
ISTART  and  NAME  vectors  are  the  same  for  the  A  and  B  matrices  and  hence 
the  comparisons  given  in  Table  5  are  conservative  for  systems  with  more 
than  one  coefficient  matrix. 

The  algorithm  to  assemble  the  element  matrices  into  the  CC  vector 
is  straight  forward  and  can  be  accomplished  in  a  dozen  or  so  steps.  In 
our  work  it  was  found  useful  to  construct  and  have  the  computer  punch 
out  the  ISTART  and  NAME  vectors  from  the  element-system  connectivity  of 
the  discretized  model  under  consideration.  Again  this  pre-processing 
algorithm  was  quite  simple  to  program. 

Solution  Methods 

It  was  not  the  intention  here  to  compare  the  many  techniques  avail- 
able for  the  solution  of  Eq.   (2),  but  rather  to  establish  the  computational 
features  of  the  multi-step,  predictor-corrector  Gear  method  for  stiff 
systems,  and  the  single-step,  implicit  Crank-Nicolson  method.     The  work 
led  to  the  development  of  a  modified  implicit  form  of  Gears  method  which 
incorporated  some  of  the  attractive  features  of  each  system. 

Having  been  specifically  developed  for  stiff  systems,  Gears  method 
handled  Eq.    (2)   quite  efficiently  with  respect  to  CPU  time.     Indeed,  Gears 
method  achieved  a  specified  accuracy  criterion  with  far  fewer  integration 
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steps  than  the  Crank-Nicolson  method,  which  experienced  particular  diffi- 
culty during  the  early  transient  period  of  time.     However,  Crank-Nicolson 
performed  equitably  with  Gears  method  during  the  later  transient  period 
as  shown  in  Figure  2. 

A  severe  drawback  to  multi-step  methods  such  as  Gears  method,  in 
contrast  to  single-step  methods  such  as  Crank-Nicolson, is  the  need  to 
retain  previous  time  solutions.     The  predictor-corrector  Gear  method,  as 
typified  in  the  DVOGER  subroutine  (5)   is  capable  of  6       order  interpolation 
in  time  and  requires  N(N  +  17)  words  of  work  storage,   in  addition  to  8N 
words  for  storage  of  past  history,  unquestionably  a  considerable  expense. 

Another  disadvantage  of  Gears  method  is  the  need  to  transform  the 
implicit  system  of  Eq.  (2)  into  the  explicit  system 

i  =  A-1    (B  -  C^  =  FM  (8) 

in  the  case  of  direct  treatment  of  the  nonlinear  term,  or 

i   =  A"1  (B  -  C*)i>  (9) 

when  the  nonlinear  term  is  given  the  linear  approximation  of  Eq.  (3).  Equa- 
tions (8)  and  (9)  are  formal  statements  since  actual  pre-multi plication 
by  A"  is  not  necessary.  These  equations  may  be  obtained  by  a  Cholesky 


decomposition  into 


LLT<1-  =  B^  -  Ci>2  (10) 


or, 

LLT4>  =  (B  -  C)i|;  (11) 

which,  in  turn,  may  be  solved  by  forward  and  backward  substitutions  for 
i.     It  should  be  noted  that  Cholesky  decomposition  eliminates  the  possible 
use  of  OCS,  and  hence  a  severe  penalty  in  coefficient  array  storage  is 
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incurred  for  the  A  matrix.  This,  in  addition  to  the  need  to  store  the 

3  p-i 

NxN  Jacobian  matrix  -r—  gives  a  great  advantage  to  the  Crank-Nicolson 
method,  despite  the  superior  ability  of  Gears  method  for  stiff  systems. 
The  modified  Gears  method,  to  be  described  shortly,  offers  a  substantial 
diminishment  in  the  storage  requirement  of  Gears  method,  while  maintain- 
ing the  advantage  of  superior  integration  for  stiff  systems.  The  net 
effect  is  a  method  which  compares  favorably  with  the  Crank-Nicolson 
method. 

The  Modified  Gear  Method 

Our  initial  work  (]_)  presented  a  dilemna  which  required  additional 
investigation.  On  the  one  hand  the  superior  tracking  ability  of  Gears 
method  for  the  stiff  system  of  equations  (2)  was  admired,  while  on  the 
other  hand  there  was  dissatisfaction  with  the  large  core  storage  require- 
ment. One  way  to  minimize  this  difficulty  was  a  modification  of  Gears 
method. 

1.  To  treat  the  implicit  system  of  Eq.  (2)  and  thus  achieve  the 
core  saving  advantages  of  OCS,  and 

2.  Eliminate  the  need  for  storing  the  NxN  Jacobian  matrix. 

In  order  to  obtain  an  understanding  of  the  modified  Gear  method 
the  following  presentation  is  given. 

The  core  storage  requirements  of  the  multi-step  Gear  method  arise 
primarily  from  the  necessity  to  compute  and  store  the  NxN  Jacobian  matrix 
3Fi/a^ . .  While  A,  B  and  C  (or  C*)  are  sparse  banded  matrices,  the 

Jacobian  A   (B-2C^),  or  A  (B-C*),  is  full.  Storage  of  this  matrix 

2 
accounts  for  N  of  the  work  storage  words  required  by  DVOGER.  At  the 
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suggestion  of  Prof.  C.  W.  Gear  (6),  an  existing  program  (7J  for  implicitly 
defined  stiff  differential  equations  was  modified  to  take  maximum  advantage 
of  the  sparse  nature  of  matrices  A,  B,  and  C  (or  C*) .  This  includes  the 
use  of  OCS  for  the  equivalent  of  the  Jacobian  matrix,  which  uses  the  same 
ISTART  and  NAME  vectors  as  the  A  and  B  matrices,  and  the  use  of  iteration 
(SOR)  to  solve  for  the  difference  between  the  predicted  value  and  the 
corrected  value  at  each  time  increment. 

A  brief  background  concerning  Gears  method  for  implicit  stiff  equations 
is  necessary  in  order  to  understand  why  SOR  results  in  a  substantial  increase 
in  computational  efficiency  for  problems  of  this  type.  In  the  general  case, 
one  seeks  to  obtain  a  numerical  solution  of  the  system  of  equations 

F  (t,  y,  y)  =  0  (12) 

with  initial  conditions  y(t  )  specified.  At  each  time  step  here  denoted 
by  superscript  k,  the  corrector  equations 

f  (t,  y(k),  -v(k)>V  +  Ek}  =  °        (13) 

(k) 
which  result  from  approximating  y  by  a  linear  combination  of  yv  '   and 

(k) 
previous  solution  values,  must  be  solved  for  yv  .  In  the  above,  h  is  the 

stepsize,  a     and  a    depend  on  the  order  of  the  method  (interpolation  on 

time),  and  z.    depends  on  the  order  and  previous  solution  values.  Gears 

Ik) 
method  uses  Newtons  method  to  solve  for  yv  ,  hence  the  Jacobian, 

J  =  - —  ("snJ  t~  >  (°r  a  reasonable  approximation  to  it)  must  be  computed. 

°     y  ik) 

Because  the  predictor  provides  a  good  estimate  of  yv    '  and  Newton's  method 
solves  for  the  difference  in  iterates,  Sy,  the  linear  equations  to  be 
solved  are  of  the  form 

J(5y)   =  -F.  (14) 
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When  the  solution  6y   is  added  to  the  predicted  value  it  is  apparent  that 
&y   need  be  accurate  to  only  a  few  significant  digits.  Also,  because  up 
to  three  Newton  iterations  are  performed,  inaccuracies  in  the  first 
iteration,  where  the  largest  corrections,  6y,   are  expected,  are  auto- 
matically accounted  for  in  succeeding  iterations.  It  should  be  noted 
that  the  iteration  used  by  Gear  is  actually  a  modi  fed  Newton  method  since 
J  is  not  evaluated  at  each  iteration,  nor  even  at  each  time  increment. 

The  error  tolerance  in  solving  for  Sy  by  SOR  was  satisfied  when 
either  of  two  criteria  were  met 

1.  The  difference  in  components  of  successive  iterates  were  as  small 

compared  to  the  solution  6y   as  that  allowed  by  Gears  method  for 

(k) 
the  relative  error  in  yv  ',  or 


2.  The  difference  in  components  of  successive  iterates  compared 

•e  1 
,(k) 


(k) 
to  yv  '   were  less  than  1  percent  of  the  error  allowed  by  Gears 


method  in  yy 

(k) 
Considering  the  comparitive  magnitudes  of  the  vectors  yv  '  and  3y,  the 

result  was  that  yery   few  SOR  iterations  are  required  to  solve  equation 

(14).  In  the  problem  under  investigation  here,  3-4  iterations  were 

required  for  the  first  Newton  iterate,  and  1-2  for  second  and  (if  required) 

third  Newton  iterates. 

A  comparison  of  the  number  of  arithmetic  operations  required  for 
solution  of  the  corrector  equation  (13)  for  three  different  schemes 
follows.  An  operation  is  defined  to  be  a  multiplication  or  division 
and  an  addition  or  subtraction. 

1.  Use  of  SOR  in  the  modified  Gear  method  results  in  M  +  6N  opera- 
tions per  iteration  of  SOR,  including  those  required  for 
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convergence  tests.     An  overhead  of  N  divisions  for  the  initial 
estimate  of  the  solution  (taken  to  be  Sy./vJ  —  )   is  required  for 
each  Newton  iteration.     Storage  for  J  requires  only  M  locations 
for  the  coefficients  since  the  ISTART  and  NAME  vectors  are  the 
same  as  for  the  A  and  B  matrices. 

2.  Use  of  Cholesky  decomposition  to  factor  J  and  solution  of  equa- 
tions (14)  by  forward  and  backward  substitutions  requires 

2Nn  -  n(n-l)  operations  for  each  Newton  iteration,  plus  an 

2 
overhead  of  about  Nn  operations  for  Cholesky  decomposition  each 

time  J  is  computed.  While  J  is  only  recomputed  eyery   few  time- 
steps,  for  large  problems  the  decomposition  of  J  will  probably 
take  as  many  operations  as  all  the  succeeding  solutions  of 
equation  (14).  The  matrix  J  is  assumed  to  be  stored  in  symmetric 
band  storage  form,  requiring  Nn  storage  locations. 

3.  In  a  routine  such  as  DVOGER  where  the  explicit  differential 
equations  must  be  used,  the  banded  nature  of  the  Jacobian  is 
destroyed  and  the  resulting  J  matrix  is  full  and  nonsymmetric. 
Gauss  elimination  for  the  solution  of  (14)  requires  N  operations 

for  forward  and  backward  substitutions,  with  the  LU  decomposition 

3 

requiring  about  N  /3  operations  each  time  J  is  computed.  Since 

2 
J  is  full,  N  storage  locations  are  required  for  it. 

The  total  amount  of  work  involved  in  solving  the  corrector  equation 
(13)  is  difficult  to  compare  in  the  first  two  instances  because  of  a 
varying  number  of  iterations  required  in  solving  (14)  by  SOR,  as  well  as 
the  variable  nature  of  M  compared  to  Nn.  It  is  clear  that  DVOGER  and 
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similar  codes  which  must  have  the  explicit  equation,  i.e.,  (8)  or  (9), 
to  solve,  cannot  compare  favorably  with  the  modified  Gear  method  which 
attacks  the  implicit  equation  (2),  since  both  storage  and  run  time  re- 
quirements are  greater.  For  the  problems  under  consideration  here,  typ- 
ical test  runs  indicated  the  solution  of  (14)  by  SOR  is  much  faster  than 
elimination,  since  the  total  run  time  for  solution  of  the  differential 
equation  was  less  by  a  factor  of  about  2.5.  With  regard  to  storage,  the 
advantage  of  SOR  over  Cholesky  decomposition  is  greater  than  shown  in 
Table  5  since  the  ISTART  and  NAME  vectors  are  the  same  as  for  the  A  and 
B  matrices.  The  details  are  presented  in  (8). 
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